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Abstract. 

The nuclear force acting between protons and neutrons is studied in the Monte 
Carlo simulations of the fundamental theory of the strong interaction, the quantum 
chromodynamics defined on the hypercubic space-time lattice. After a brief summary 
of the empirical nucleon-nucleon (NN) potentials which can fit the NN scattering 
experiments in high precision, we outline the basic formulation to derive the potential 
between the extended objects such as the nucleons composed of quarks. The equal- 
time Bethe-Salpeter amplitude is a key ingredient for defining the NN potential on the 
lattice. We show the results of the numerical simulations on a 32^ lattice with the 
lattice spacing a ~ 0.137fm (lattice volume (4.4 fm)"*) in the quenched approximation. 
The calculation was carried out using the massively parallel computer Blue Gene/L 
at KEK. We found that the calculated NN potential at low energy has basic features 
expected from the empirical NN potentials; attraction at long and medium distances 
and the repulsive core at short distance. Various future directions along this line of 
research are also summarized. 



PACS numbers: 12.38.Gc, 13.75.Cs, 21.30.-x 



Nuclear Force from Lattice QCD 



2 



1. Introduction 

One of the long standing problems in particle and nuclear physics is the origin of the 
strong nuclear force which holds the nucleons (protons and neutrons) inside atomic 
nuclei. For the past half century, phenomenological fits of the proton-proton (pp) 
and neutron-proton (np) scattering data assuming empirical nucleon-nucleon (NN) 
potentials have been attempted [1, 2]: The potentials, which can fit more than 2000 
data points of the NN phase shift with xV^of ~ 1 for Ti^b < SOOMeV, include the CD- 
Bonn potential [3], Argonne Vis potential [4] and Nijmegen potentials [5]. Alternative 
approach on the basis of the chiral perturbation theory has been also developed [6]. 

Shown in Fig.l are three examples of the empirical NN potentials in the ^5*0 
channel.^ From this figure, some characteristic features of the nuclear force can be 
seen: 

I. The long range part of the nuclear force (r > 2 fm) is dominated by the one pion 
exchange originally introduced by Yukawa [8]. 

II. The medium range part (1 fm < r < 2 fm) receives significant contributions from 
the exchange of multi-pions and heavy mesons (p, cu, and a). In particular, the 
spin-isospin independent attraction of about —50 ~ —100 MeV in this region plays 
an essential role for binding the atomic nuclei. 

III. The short range part (r < 1 fm) behaves as a repulsive core originally introduced by 
Jastrow [9] to explain the pp and np scattering phase shifts simultaneously. Such a 
short range repulsion is relevant to the stability of atomic nuclei, to the maximum 
mass of neutron stars, and to the ignition of the Type II supernova explosions [10]. 

It is now well established that the nucleons are made of quarks and gluons which 
obey the law of quantum chromodynamcs (QCD) [11]. Therefore, it is tempting to 
derive the strong nuclear force from the quark-gluon degrees of freedom. Two nucleons 
have a sizable overlap at distance r < 1 fm since the proton's charge radius is about 
0.86 fm. Therefore, the nuclear force at short distances must be described by taking 
into account the direct exchange of quarks and gluons between the nucleons. So far, 
there have been numerous theoretical attempts to understand the nuclear force from the 
quark structure of the nucleon [12]. However, conclusive results have not been obtained 
because of the highly non-perturbative nature of QCD. 

Recently, the present authors have developed a new approach to this long standing 
problem of nuclear force on the basis of QCD defined on a space-time lattice (LQCD) 
[13, 14]. In LQCD, physical quantities are expressed by highly multi-dimensional 
integrals which can be carried out by importance sampling method. Our idea is to 
define the NN potential in the coordinate space from the equal-time Bcthe-Salpeter 
amplitude calculated on the lattice: The first results on the central NN potential in the 
^5*0 and ^5"! channels have been reported with the quenched lattice QCD simulations in 

I A system of two nucleons with total spin s = (0, 1), orbital angular momentum L = {S, P,D,- ■ •) 
and total angular momentum J = (0, 1, 2, • • •) is denoted as 
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Figure 1. Three examples of the high-precision NN potentials in the So channel. 
AV18 stands for the Argonne uig potential [4], Reid93 stands for one of the Nijmegen 
potentials [5] and Bonn stands for the Bonn potential [7]. I, II and III correspond 
to the long range part, medium range part and the short range part, respectively, as 
discussed in the text. 

[15, 16, 17]. Also the first results of the hyperon-nucleon potential have been reported 
in [18, 19]. 

In the following, we will outline the field theoretical derivation of the NN potential 
from QCD [20] in Sec. 2 and Sec. 3. Then, we show how to define the potential in LQCD 
formalism in Sec. 4. The basic setup and the method of our numerical simulations are 
shown in Sec. 5. Some numerical results of the low energy NN potential taken from 
[15, 16] are shown in Sec. 6. The last section is devoted to summary and concluding 
remarks with a discussion on the future directions. 

2. Bethe-Salpeter wave function and the NN potential 

Let us start with the definition of the Bethe-Salpeter (BS) amplitude for the proton- 
neutron system, 

'^apix.y) = {0\T[np{y)pa{x)]\p{q,s)n{q',s');m), (1) 
fif^iy) = eabc{ua{y)C-f5db{y)) dc/siy), (2) 

Pa{x) = eabc{ua{x)C-f5db{x)^Uca{x), (3) 

where (g, s) and (g ', s') denote the spatial momentum and the spin-state of the incoming 
proton and those of the neutron, respectively. The local composite operator for the 
neutron (proton) are denoted by fip^y) {pa{x)) with the operators for the up-quark u{x) 
and the down-quark d{x). Also, a and /5 denote the Dirac indices, a, b and c the color 
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indices, and C the charge conjugation matrix in the spinor space. 

One of the advantages to use local operator for the nucleon is that the Nishijima 
and Zimmermann (NZ) reduction formula [21] for local composite fields can be utilized. 
In particular, in and out nucleon fields are defined through the Yang-Feldman equation, 

VZNm/out{x) = N{x) - j 5'adv/ret(2^ " x' ] m^) J {x')dx' , (4) 

where N{x) is the nucleon composite operator {n{x) or p{x) in Eqs.(2,3)), iVin/out(a;) is 
the associated in/out field, mN is the physical nucleon mass, and J{x) = {i(^x—T^N)N{x). 
The advanced/retarded propagator for the free nucleon field with the mass mN is denoted 
by Sady/ret{x — x';mN). The normalization factor, y/Z, is the coupling strength of the 
composite operator N{x) to the physical nucleon state. 

Through the NZ reduction formula, the BS amplitude in Eq.(l) is related to the 
four-point Green's function G4 of the composite nucleons which is decomposed into the 
free part and the scattering part, G4 = Z'^{Gf^ +04^^^), where Cf^ is proportional to a 
product of the free nucleon propagators. After taking the equal time limit, x^ = = t, 
with f = X — y, it is straight forward to rewrite Eq.(l) to the the following integral 
equation in the cm. frame (g ' = —q) [20], 



"^a^ir, t) = ^^^{f- q, s, s')e-'^V?^\ (5) 
i>a/3{r; s, s') = Zu^i^q, s)up{-q, s')e"^"' 

7^—^ e''''^J^a/3;js{k; q)u^{q, s)us{-q, s'). (6) 



Here '^^(g; s) is the positive-energy plain- wave solution of the Dirac equation and 
^ais-^sik] (f) is an integral kernel obtained from G^^^ after carrying out the /cQ-integration. 
Hereafter, we call ■ipa/di'f^lQ, s, s') as the Bethe-Salpeter wave function. The differential 
form of the above equation is obtained by multiplying (g ^ + V^)/mN to Eq.(6); 

— {q^ + ^^)ipapif; q, s, s') = K^p{r, q, s, s'). (7) 
mN 

Note that we have not made any non-relativistic approximation to derive Eq.(7). An 
important observation is that the plain wave component of ipa/si'^', Q, s, s') is projected out 
by the operator (g ^ + V^) so that the function Kap{r; q, s, s') is localized in coordinate 
space as long as |gj stays below the inelastic threshold as noted for pion-pion scattering 
in [22]. This is equivalently said that the Fourier transform of K with respect to r, 
which is proportional to the half off-shell T-matrix relating the the on-shell state with 
momentum q and the off-shell state with momentum k, does not develop a real pole as 
a function of if \q\ is below the inelastic threshold. 

3. The NN potential 

In an abbreviated notation, (r, a, (3) x, and {q, s, s') q, Eq.(7) is written as 

{E,-HoMx;q)^K{x;q), (8) 
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where Eg — (f /m^ and Hq = — V^/mN- This equation defined in a finite box can be used 
to extract various information on the NN scattering from the lattice QCD simulations: 

(i) Consider K{x; g) as a measure to identify the length R beyond which the two 
nucleons do not interact. If we stay in such a region where K{x > R]q) — 0, 
the wave function ip{x > R; q) can be expanded by the solution of the Helmholtz 
equation inside a finite box. Then one can extract the phase shift given the incoming 
energy Eg. This is the approach originally proposed in [23] and is later elaborated 
to study hadron-hadron scatterings on the lattice [24, 22, 25, 26]. 

(ii) Alternatively, one may extract the half off-shell T-matrix in momentum space by 
calculating the left-hand-side of Eq.(8) in the coordinate space and making Fourier 
transform with respect to x. 

(iii) One can go one-step further and define the non-local NN potential U (x, x') from 
K{x; g), so that Eq.(8) becomes the Schrodinger type equation. 

If we are interested only in the NN scattering phase shift in the free space, the 
procedure (i) is certainly enough. On the other hand, if we are interested in applying 
Eq.(8) to the problems of bound states and the nuclear many-body system, (ii) and 
(iii) are useful since they give us the off-shell information in a well-defined manner in 
QCD. To see this explicitly for the case (iii), we introduce a set of functions labeled by 
q, {ip{x; q)}, which is dual to the set {ip{x; q)} in the following sense: 

J dx tpix; q)i/j{x;p) = 5g^p. (9) 

As long as the dimensions of the x-space and p-spacc arc the same and the elements 
in {ip{x;p)} are linearly independent, such a dual basis exists and is unique. If the 
dimension of p-space is less than that of x-space, the dual basis exists but is not unique. § 
Assuming the existence (but not necessarily the uniqueness) of the dual basis, the non- 
local potential can be defined as 

U{x,x') = J dp K{x,p)ijj{x;p). (10) 

The Eqs.(9,10) lead to the formula, K{x;q) — J dx' U{x,x')il}{x';q), so that Eq.(8) 
becomes 

Va/3(r; q, s, s') + / d^r' Uap;^s{r, f ')ip^s{r '; q, s, s') = Egijjap{f : q, s, s'). (11) 

77T.N J 

Note that the non-local potential U can be rewritten in the form, 

U{f,f') = V{f,V)6{f-f'). (12) 

§ This is easily seen as follows. Let us introduce a basis {e(i), 6(2), • • • , e(^')} in the A''-dimensional 
vector space. The BS wave function in discretized coordinates ip{x;p) = ■(/'(*; corresponds to e^^^ 
with 1 <i < N and 1 < a < M < N. If M = N, there exits a unique dual basis, {e(i),e(2), • • • ,e(jv)} 
satisfying e(Q,) • 6(^3) = 6af3 (see any textbook of linear algebra). If M < A^, there is still a dual basis 
{e(i),e(2), • • • ,e(M)} satisfying the above condition for a < M and (3 < M. However, it is not unique 
because one always has a freedom to add linear combinations of e-^ (M + 1 < 7 < A'') to the above dual 
basis. 
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The general structure of V{f, V) under various symmetry constraints in the non- 
relativistic kinematics has been worked out by Okubo and Marshak [27]. If we further 
make the derivative expansion at low energies [28] , we obtain the expression familiar in 
the phenomenological potentials acting on the upper components of the wave function; 

V{f, V) = Voir) + K(r)(ai ■ a^) + K(r)(ri ■ fs) + Kr(r)(ai ■ a2)(ri ■ T2) 

+ VT{r)S,2 + VTr{r)S,2{ri ■ fs) + Vis(r)(L • S) + Vi^sr{r){L ■ S){n ■ T2) 
+ 0(V2). (13) 

— * 

Here S12 — 3{ai-n){a2-n) —Bi-a2 is the tensor operator with ft — r/\r\, S — (cti + CT2)/2 

— * 

the total spin operators, L = —ifx V the orbital angular momentum operator, and ri 2 
are the isospin operators for the nucleons. Each component of the potential in Eq.(13) 
can be obtained by appropriate spin, isospin and angular momentum projection of the 
BS wave function. Also, the higher derivative terms of the potential in Eq.(13) can be 
deduced by combining the BS wave functions for different incident energies. 

It is in order here to remark that the structure of the non-local potential U {x, x') is 
directly related to the nucleon interpolating operator adopted in defining the Bethe- 
Salpeter wave function. Different choices of the interpolating operator would give 
different forms of the NN potential at short distance, although they give the same 
phase shift at asymptotic large distance. The advantage of working in QCD is that we 
can unambiguously trace the connection between the NN potential and the interpolating 
operator. 



4. Effective central potential on the lattice 

The BS wave function in the S-wave on the lattice with the lattice spacing a and the 
spatial lattice volume is obtained by 

^(r) = ^ E E ^^/^ (0 \Mm + ^)Paix)\pn; q) . (14) 

The summation over i? G O is taken for the cubic transformation group to project out 
the iS-wavc. II The summation over x is to select the state with zero total momentum. We 
take the upper components of the Dirac indices to construct the spin singlet (triplet) 
channel by p^-^™siet _ (^g-^^^^ (-Pa/3~*"'''^* ~ (o"i)a/3)- The BS wave function tjj{f) is 
understood as the probability amplitude to find "neutron-like" three-quarks located at 
point X + f and "proton-like" three-quarks located at point x. 

In the actual simulations, the BS wave function is obtained from the four-point 
correlator, 

G4{x,y,t;to) = (o |n;3(y, i)Pa(^, *) Jpn(to)| o) = A (0 |n^(y)p«(^)| ri) e-^"(*-*°). (15) 

n 

Here j7'p„(to) is a wall source located at t = to; which is defined by J'pn{to) = 
Pap Y^x,yPa{x, to)hi3{y, to). The eigen-energy and the eigen-state of the six quark system 

II In principle, this projection cannot remove possible contamination from the higher orbital waves with 
L > 4, although these contributions are expected to be negligible. 
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are denoted by En and |n), respectively, with the matrix element A„(to) = {n\J^pn{to)\0). 
For {t — to)/a ^ 1, the G4 and hence the wave function ip are dominated by the lowest 
energy state. 

The lowest energy state created by the wall source j7'p„(to) contains not only the 
5'-wave component but also the D-wave component induced by the tensor force. In 
principle, they can be disentangled by preparing appropriate operator sets for the sink. 
Study along this line to extract the mixing between the 5'-wave and the D-wave at 
low energies has been put forward recently in [29]. In the present paper, instead of 
making such decomposition, we define an "effective" central potential Vc(r) according 
to Refs.[15, 16]: 

Vc{r)^E, + —-^^-p^. (16) 

Note that one can test the non-locality of the potential U {x, x') by evaluating the 
effective central potential for different energies. If there arises appreciable energy 
dependence in Vc(r), it is a signature of the necessity of high derivative terms in Eq.(13). 



5. Setup of the lattice simulations 

In lattice QCD simulations, the vacuum expectation value of an operator 0{q, q, U) is 
defined as 

(O) = / n dU{€) n dq{x)dq{x) 0{q, q, U) e-5ffe?"'f^)-5s(^) (17) 
= Z-^ fl[dU{i) Q{U)detM{U) e'^^^^l (18) 

where Z = J UedU{e)Uxdq{x)dq{x) e^^f^i^^^^'^-^eiu) jg ^^e QCD partition function, 
Sf = J2x,x' ^xx' {U)q{x') is the quark part of the action, and Sg is the gluon part of 
the action. The quark field q{x) is defined on each site x of the hypercubic space-time 
lattice, while the gluon field U (i) denoted by 3 x 3 special unitary matrix is defined on 
each link i. In Eq.(18), the integration of the quark fields is carried out analytically. In 
the quenched approximation adopted in our simulation, the virtual fermion loop denoted 
by det M{U) is set to be 1 and the integration over link variables U is performed using 
the importance sampling method [30]. 

We employ the standard plaquette gauge action on a 32*^ lattice with the bare QCD 
coupling constant P = 6/ = 5.7. The corresponding lattice spacing is determined to be 
1/a — 1.44(2) GcV [a ~ 0.137 fm) from the p meson mass in the chiral limit [31]. Then, 
the physical size of our lattice becomes L ~ 4.4 fm. For the fermion action, we adopt 
the standard Wilson quark action with the hopping parameter (k, — 0.1640, 0.1665 and 
0.1678) which controls the quark masses. The periodic (Dirichlet) boundary condition 
is imposed on the quark fields along the spatial (temporal) direction. To generate the 
quenched gauge configurations, we adopt the heatbath algorithm with overrelaxation 
and sample configurations are taken in every 200 sweeps after skipping 3000 sweeps for 
thermalization. 
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K 


Nconf 


[MeV] 


mN [MeV] 


(t - to)/a 


Egi^So) [MeV] 


E.i^Si) [MeV] 


0.1640 


1000 


732.1(4) 


1558.4(63) 


7 


-0.400(83) 


-0.480(97) 


0.1665 


2000 


529.0(4) 


1333.8(82) 


6 


-0.509(94) 


-0.560(114) 


0.1678 


2021 


379.7(9) 


1196.6(32) 


5 


-0.675(264) 


-0.968(374) 



Table 1. The number of gauge configurations A^conf, the pion mass rriTr, the nucleon 
mass rriN, time-shce t — to on which BS wave functions are measured, and the non- 
relativistic energies Eg = q^/m^ for ^Sq and channels. 



For our numerical simulations, we use IBM Blue Gene/L at KEK, which consists of 
10,240 computation nodes with total theoretical performance of 57.3 TFlops. A modified 
version of the CPS++ (the Columbia Physics System) [33] is used to generate quenched 
gauge configurations and propagators of quarks. Most of the computational time is 
devoted to the calculation of the four-point function of nucleons, for which our code 
achieves 34-48 % of peak performance. Totally about 4000 hours are used by queues 
with 512 nodes for the calculations of the effective central potentials corresponding to 
the three values of hopping parameters. The number of sampled gauge configurations 
A^conf, the pion mass tHt^, and the nucleon mass ttin are summarized in Table 1. (For 
K — 0.1678, we have removed 24 exceptional gauge configurations from the sample.) 

We adopt the wall source on the time-slice t/a = to/a = 5. The BS wave functions 
are measured on the time-slice {t — to)/a = 7,6,5 for k, = 0.1640,0.1646,0.1678, 
respectively. The ground state saturation is examined by the t-dependence of the NN 
potential. We employ the nearest neighbor representation of the discretized Laplacian 
as V^/(a;) = J2f^i {f{x + afii) + f{x — an,)} — 6f{x), where rij denotes the unit vector 
along the i-th coordinate axis. BS wave functions are fully measured for r < 0.7 fm, 
where rapid change of the NN potential is expected. Since the change is rather modest 
for r > 0.7 fm, the measurement of BS wave functions is restricted on the coordinate 
axes and their nearest neighbors to reduce the calculational cost. The "asymptotic 
momentum" q is obtained by fitting the BS wave function with the Green's function in 
a finite and periodic box [23]: 

1 i(27r/L)n-r 

which satisfies (V^ + q'^)G{f; q^) = —SL{r) with Silr) being the periodic delta- function. 
In the actual calculation, Eq.(19) is rewritten in terms of the heat kernel Ti satisfying the 
heat equation, dt7i{t,r) = V^7Y(t, r) with the initial condition, 7Y(t 0+,r) = Si^r)- 
The fits are performed outside the range of NN interaction determined by V'^ip{r) /ip{r) 
[25]. 

6. Numerical results 

Fig. 2 (upper panel) shows the BS wave functions in ^5*0 and ^5*1 channels for k = 0.1665. 
The suppression of the wave function in the region r < 0.5 fm indicate the existence of 
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Figure 2. (Upper panel) The NN wave functions in ^Sq and "^^i channels. The inset 
is a 3D plot of the wave function (f){x, y,z = 0; ^Sq). (Lower panel) The NN effective 
central potential in the ^Sq and ^5*1 channels for m,r = 529 MeV (k — 0.1665). 



repulsion at short distance. Fig. 2 (lower panel) shows the reconstructed NN potentials 
for K = 0.1665, i.e., the effective central potentials for ^Sq and channels. See Table 1, 
for the values of the non-relativistic energies Eg = q'^/m^ in Eq.(16). 

We show the NN potentials for three different quark masses in ^5*0 channel in Fig. 
3. As the quark mass decreases, the repulsive core at short distance is enhanced rapidly, 
whereas the attraction at medium distance is modestly enhanced. This indicates that it 
is important to perform the lattice QCD calculations for the lighter quark mass region 
in order to make quantitative comparison of our results with the observables in the real 
world. 

Although there exist both attraction and repulsion, the net effect of our potential 
is attractive at low energies as shown by the scattering lengths calculated from 
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Figure 3. Central potentials in the 5*0 channel for three different quark masses in 
the quenched QCD simulations on a (4.4fm)'' lattice. 





[MeV] 


00(^50) [fm] 


aoCSi) [fm] 


0.1640 


732.1(4) 


0.115(26) 


0.140(31) 


0.1665 


529.0(4) 


0.126(25) 


0.140(31) 


0.1678 


379.7(9) 


0.159(66) 


0.252(104) 



Table 2. Scattering lengths obtained from the Liischer's formula [23] in the spin- 
singlet and spin-triplet channels for different quark masses. 



the Liischer's formula in Table 2 [23]. The attractive nature of our potential is 
qualitatively understood by the Born approximation formula for the scattering length 
Co — —rriM J Vc{r)r'^dr. Owing to the volume factor r'^dr, the attraction at medium 
distance can overcome the repulsive core at short distance. However, there is a 
considerable discrepancy between the scattering lengths in Table 2 and the empirical 
values, ao'^''^^(^5'o) ~ 20 fm and ao'^''^^(^5'i) ~ —5 fm. This is attributed to the heavy 
quark masses employed in our simulations. 

If we can get closer to the physical quark mass, there appears an "unitary region" 
where the NN scattering length becomes singular and changes sign as a function of 
the quark mass [31, 32, 26]. The singularity is associated with the formation of the 
di-nucleon bound state. Because of this, the NN scattering length becomes a highly 
non-linear function of the light quark mass. One should note that the NN potential 
changes smoothly even in the unitary region in contrast to the scattering length. This 
is one of the reasons why the NN potential is much more appropriate quantity to be 
examined on the lattice instead of the NN scattering length. 
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7. Summary and concluding remcirks 

In this article, we have outhned the basic notion of the nucleon-nucleon potential and 
its field-theoretical derivation from the equal-time Bethe-Salpeter amplitude. Such a 
formulation allows us to extract the potential between extended objects by using the 
lattice QCD simulations. The central part of the NN potential at low energies was 
obtained in lattice QCD simulations with quenched approximation. It was found that 
the NN potential calculated on the lattice at low energy shows all the basic features 
expected from the empirical NN potentials determined from the NN scattering data; 
attraction at long and medium distances and the repulsion at short distance. This is 
the first step toward the understanding of atomic nuclei from the fundamental law of 
the strong interaction, the quantum chromodynamics. 

There are a number of directions to be explored on the basis of our approach: 

1. Energy dependence of the NN potential in Eq.(16) should be studied to test the 
non-locality of the potential U {x, x') and the validity of its derivative expansion. 
This is currently under investigation by changing the spatial boundary condition 
of the fermion field [34] . 

2. The tensor force, which mixes the states with different orbital angular momentum 
by two units, is a unique feature of the nuclear force and plays an essential role 
for the deuteron binding. This is also under investigation by projecting out the 

^Si component and ^Di component separately from the exact two-nucleon wave 
function with J = 1 on the lattice [29, 35]. The spin-orbit force, which is known to 
be strong at short distances in empirical NN force, should be also studied. 

3. Three nucleon force is thought to play important roles in nuclear structure and 
also in the equation of state of high density matter. Since the experimental 
information is scarce, simulations of the three nucleons on the lattice combined 
with appropriate generalizations of the formulas in Section 2 may lead to the 
first principle determination of the three nucleon potential in the future. With 
these generalizations of the present approach, one may eventually make a firm link 
between QCD and the physics of nuclear structure [36]. 

4. The hyperon-nucleon (YN) and hyperon-hyperon (YY) potentials are essential for 
understanding the properties of hyper nuclei and the hyperonic matter inside the 
neutron stars. However, the experimental data are very limited due to the short 
life-time of hyperons. On the lattice, NN, YN and YY interactions can be treated 
in the same footing since the difference is only the mass of the strange quark. 
Recently, the SN potential [18, 19] and the AA^ potential [37] are examined as a 
first step toward systematic derivation of the hyperon interactions. 

5. To compare the NN potential on the lattice with experimental observables, it 
is necessary to carry out full QCD simulations which take into account the 
dynamical quark loops. Study of the nuclear force with the use of the fuU- 
QCD configurations generated by PACS-CS Collaboration [38] is currently under 
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investigation [34, 35, 37]. 
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